/* Declaration of set of functions to solve numerictly */
#include "yyt05Set.h"
using namespace std;

/* FROM YYT - GA METOD */

yyt::yyt(Parameters *Par, int numberOfEquation)
{
	// cerr << "konstruktor yyt: this " << this << endl;
	par = Par;
	sdelta = (*Par)[0]; // 0. sdelta = 1-n (small delta),
	J = (*Par)[1]; // 1. J (exchange coupling constant)
	t = (*Par)[2]; // 2. t (hopping integral for the nearest neighbors),
	U = (*Par)[3]; // 3. U (repulsive potential for electrons in the same spatial point).
	beta = (*Par)[4]; // 4. beta (inverse temperature),
	L = (*Par)[5]; // 5. L (siza of lattice)
	nEqu = numberOfEquation;
	yytF = new yytFunctions(Par);
}

yyt::~yyt()
{
	// cerr << "destruktor yyt" << endl;
	delete yytF;
}

yytGA::yytGA(Parameters *Par, int _nEqu, Parameters *X0) : yyt(Par, _nEqu)
{
	// cerr << "konstruktor yytGA: this " << this << endl;
	setX0(X0);
	// cerr << "DEBUG in yytGA: " << nEqu << endl;
}

void yytGA::setX0(Parameters *X0)
{
	this->X0 = X0;
	yytF->setX0(X0);

	gs = yytF->gs();
	dgtdm = yytF->dgtdm();
	dgtdd = yytF->dgtdd();
	dgsdm = yytF->dgsdm();
	dgsdd = yytF->dgsdd();
	BDelta_d = yytF->BDelta_d();
	BDelta_af = yytF->BDelta_af();
}

void yytGA::onePoint(double *y, int i, int j)
{
	enum Variables { BDelta, chi, m, d, mu };

	double gamma = yytF->gamma(i*2.*M_PI/L, j*2.*M_PI/L);
	double eta = yytF->eta(i*2.*M_PI/L, j*2.*M_PI/L);
	double epsilon = yytF->epsilon(i*2.*M_PI/L, j*2.*M_PI/L);
	double dzeta = sqrt(epsilon*epsilon + BDelta_af*BDelta_af);
		double a1 = dzeta - (*X0)[mu];
		double a2 = -dzeta-(*X0)[mu];
		double b = BDelta_d*eta;
	double E1 = sqrt(a1*a1 + b*b);
	double E2 = sqrt(a2*a2 + b*b);

	double part = (dzeta-(*X0)[mu])/E1 + (dzeta+(*X0)[mu])/E2;

	y[BDelta] += eta*eta*BDelta_d*(1./E1 + 1./E2);
	y[chi] += gamma*epsilon/dzeta*((-dzeta+(*X0)[mu])/E1 + (-dzeta-(*X0)[mu])/E2);
	y[m] += BDelta_af/dzeta*part - gamma*epsilon/dzeta*part*(dgtdm*t-dgsdm*3./8.*J*(*X0)[chi])/(2.*gs*J);
	y[d] += gamma*epsilon/dzeta*part;
	y[mu] += (dzeta-(*X0)[mu])/E1 + (-dzeta-(*X0)[mu])/E2;
//	printf("<=> %20.20f %20.20f %20.20f =  %20.20f\n", yytF->gs(), J, (*X0)[BDelta], yytF->BDelta_d());
}

double* yytGA::sumationAll()
{
	enum Variables { BDelta, chi, m, d, mu };
//	sumDE();
	this->nEqu = 5;
	double *y = new double[nEqu];
	clear(y);

		double *temp = new double[nEqu]; // temporary vector using during summation
	    clear(temp);
	    unsigned int max = static_cast<int>(L)/2; // attention: L should be even!

	    // SMART WAY OF SUMMATION
	    // I use symmetry k_x,y <-> -k_x,y
	    // * edges (2 times)
	    for(unsigned int i=1; i<max; i++)
	    {
	    	onePoint(temp, i, 0);
	    	onePoint(temp, 0, i);
	    	onePoint(temp, i, max-i);
	    }
	    for(unsigned int i=0; i<nEqu; i++) { y[i] += 2.*temp[i]; temp[i] = 0; }

	    // * interior area (4 times)
	    for(unsigned int i=1; i<max; i++)
	       for(unsigned int j=1; j<i; j++)
	       {
	    	   onePoint(temp, max-i, j);
	       }
	    for(unsigned int i=0; i<nEqu; i++) { y[i] += 4.*temp[i]; temp[i] = 0; }

	    // * missing corners
	    onePoint(y, max, 0);
	    onePoint(y, 0, max);
	    onePoint(y, 0, 0);

	double chi2 = (*X0)[chi]*(*X0)[chi];
	double m2 = (*X0)[m]*(*X0)[m];
	double BDelta2 = (*X0)[BDelta]*(*X0)[BDelta];

	y[BDelta] /= 4.;
	y[chi] /= 4.;
	y[m] /= 2.;
	for(unsigned int i=0; i<nEqu; i++) { y[i] /= L*L; }
	y[BDelta] -= (*X0)[BDelta];
	y[chi] -= (*X0)[chi];
	y[m] -= (*X0)[m] + (3./4.*(BDelta2+chi2) + 2.*m2)/(2.*gs);
	y[d] *= dgtdd*t + 3./8.*J*(*X0)[chi]*dgsdd;
		y[d] += U + dgsdd*(3./4.*J*(BDelta2+chi2) + 2.*J*m2);
	y[mu] -= sdelta;

//	if((*X0)[d] < 0.) y[d] for(unsigned int i=0; i<nEqu; i++) { y[i] = 10000; }

	return y;
}
